Stability and Conditioning of a Modified Lorenz System

Shane McQuarrie
Math 510
18 December 2015

Section 1: Introduction

The Lorenz system is one of the fundamental models for chaos. We are interested in examining a modification that extends the system to five dimensions. On the left is the original system, and on the right, the modified version.

$$ \begin{align} \nonumber &\dot{X} = \sigma \left(Y-X\right) && \dot{X} = \sigma \left(Y-X\right) - pU\\ \nonumber &\dot{Y} = rX-XZ-Y && \dot{Y} = rX-XZ-Y\\ \nonumber &\dot{Z} = XY-bZ && \dot{Z} = XY-bZ\\ \nonumber & && \dot{U} = -XV + pX - \sigma U\\ \nonumber & && \dot{V} = XU - qV \end{align} $$

where $p = \frac{\sqrt{\sigma}b^{3/2}}{\epsilon\pi^2 2^{11/4}}$ and $q = (4-b)\sigma$.

We begin by importing the necessary modules.

In [1]:
%matplotlib inline
from matplotlib import rcParams, pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import odeint
from scipy.linalg import norm
import numpy as np
rcParams['figure.figsize'] = (16,10)

First we write a function for solving the original Lorenz equations. We may specifiy the system parameters $b$, $r$, and $\sigma$, as well as initial conditions, the time interval, and tolerances for the ODE solver.

For consistency, we'll use the same initial conditions repeatedly.

In [2]:
# Set initial conditions (changing this line will change the rest of the notebook).
initial = np.ones(5)*10

def solve_original(b=8/3., r=28., sig=10.,
                   x0=initial[:3], t=20, atol=1e-14, rtol=1e-12):
    """Generate a solution of the original Lorenz system."""

    # Initialize remaining parameters.
    time = np.linspace(0, t+1, t*100)

    # Define the system.
    def original(data, T):
        """Original Lorenz system."""
        X,Y,Z = data
        Xdot = sig*(Y-X)
        Ydot = r*X - X*Z - Y
        Zdot = X*Y - b*Z
        return np.array([Xdot, Ydot, Zdot])

    # Solve and return the solution.
    return odeint(original, x0, time, atol=atol, rtol=rtol).T

To visualize the solutions of the system we write a function to plot a solution set in 3 dimensions. Unfortunately, IPython Notebook doesn't cooperate with animations very well, so the images in this notebook are static.

In [3]:
def plot_single(sol, color='b', title=None):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot(sol[0], sol[1], sol[2], color=color)
    ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')
    ax.set_xlim3d([min(sol[0]), max(sol[0])])
    ax.set_ylim3d([min(sol[1]), max(sol[1])])
    ax.set_zlim3d([min(sol[2]), max(sol[2])])
    plt.title(title)

    plt.show()

plot_single(solve_original(), title="Original Lorenz System")

We similarly define a function to solve for the modified Lorenz system. Note the similarity betweeen the original and modified system solutions.

In [4]:
def solve_modified(b=8/3., r=28., sig=10., eps=.05,
                   x0=initial, t=20, atol=1e-14, rtol=1e-12):
    """Generate a solution of the modified Lorenz system."""

    # Initialize remaining parameters.
    q = (4 - b)*sig
    p = (np.sqrt(sig)*b**(1.5))/(eps*np.pi**2*2**(11./4.))
    time = np.linspace(0, t+1, t*100)

    # Define the system.
    def modified(data, T):
        """Modified Lorenz system."""
        X,Y,Z,U,V = data
        Xdot = sig*(Y - X) - p*U
        Ydot = r*X - X*Z - Y
        Zdot = X*Y - b*Z
        Udot = -1.*X*V + p*X - sig*U
        Vdot = X*U - q*V
        return np.array([Xdot, Ydot, Zdot, Udot, Vdot])

    # Solve and return the solution.
    return odeint(modified, x0, time, atol=atol, rtol=rtol).T

def plot_double(sol1, sol2, color1='b', color2='r', title=None):

    assert sol1.shape[1] == sol2.shape[1], "datasets not aligned"
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    curve1, = ax.plot(sol1[0], sol1[1], sol1[2], color=color1)
    curve2, = ax.plot(sol2[0], sol2[1], sol2[2], color=color2)
    ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')
    ax.set_xlim3d([min([min(sol1[0]), min(sol2[0])]),
                   max([max(sol1[0]), max(sol2[0])])])
    ax.set_ylim3d([min([min(sol1[1]), min(sol2[1])]),
                   max([max(sol1[1]), max(sol2[1])])])
    ax.set_zlim3d([min([min(sol1[2]), min(sol2[2])]),
                   max([max(sol1[2]), max(sol2[2])])])
    plt.title(title)
    plt.show()


s1 = solve_original()
s2 = solve_modified()

plot_single(s1, title="Original System", color='b')
plot_single(s2, title="Modified System", color='r')
plot_double(s1, s2, color1='b',color2='r',
            title="Overlay of both solutions. The Blue is the Original System and the Red is the Modified System.")

The two solutions tend to be qualitatively similar in $\mathbb{R}^3$.

Section 2: Stability

We analyze the stability of each of these systems visually (the rigorous stability of the ODE solver is a little beyond the scope of this project).

We will examine three types of perturbations: small changes in the initial condition, different tolerances in the ODE solver, and small changes in some of the system parameters.

2.1: Perturbations in the Initial Condition

Since the Lorenz system is highly chaotic, we expect that small changes in the initial conditions will result in significant changes in the solutions, provided that we allow the system to evolve for long enough.

To measure this effect, we allow the system to evolve for 3000 time intervals, and check the differences at intervals of 1000 units.

In [5]:
# Create a small perturbation and add it to the initial condition.
perturbation = np.random.randn(5)*.001
perturbed = initial + perturbation

# Solve the original and modified systems, with the original and perturbed initial conditions.
initial_org1 = solve_original(x0=initial[:3], t=30)
initial_org2 = solve_original(x0=perturbed[:3], t=30)
initial_mod1 = solve_modified(x0=initial, t=30)
initial_mod2 = solve_modified(x0=perturbed, t=30)

initial_diff1 = np.abs(np.mean(initial_org1-initial_org2, axis=0))
initial_diff2 = np.abs(np.mean(initial_mod1-initial_mod2, axis=0))

# Plot results.
plot_double(initial_org1, initial_org2, color1='b', color2='g',
            title='Original System with Perturbed Initial Conditions')
plot_double(initial_mod1, initial_mod2, color1='r', color2='g',
            title='Modified System with Perturbed Initial Conditions')

At first, the solutions appear to be identical. However, it doesn't take too long for them to separate, which is why we can see multiple colors. What's most interesting, however, is that the modified system appears to be a little less affected than the original system as time increases.

In [6]:
def analyze_difference(diff1, diff2, iters, verbose=True):
    """Print and plot statistics on the arrays `diff1` and `diff2` with increments `iters`."""
    def analyze(diff):
        """Print the maximum, minimum, and average values in an array of differences."""
        ma, me, no = np.max(diff), np.mean(diff), norm(diff)
        if verbose is True:
            print("""After {} units:
                Maximum difference:   {}
                Average difference:   {}
                2-Norm of difference: {}
            """.format(len(diff), ma, me, no))
        return ma, me, no
    
    org_maxs, org_means, org_norms = [], [], []
    mod_maxs, mod_means, mod_norms = [], [], []
    if verbose is True:
        print("Original System --------------------------------------------------")
    for stop in iters:
        x1, y1, z1 = analyze(diff1[:stop])
        org_maxs.append(x1); org_means.append(y1); org_norms.append(z1)
    if verbose is True:
        print("Modified System --------------------------------------------------")
    for stop in iters:
        x2, y2, z2 = analyze(diff2[:stop])
        mod_maxs.append(x2); mod_means.append(y2); mod_norms.append(z2)
    
    plt.subplot(121)
    plt.plot(org_maxs, linewidth=2, label="Original: Max Difference")
    plt.plot(org_means, linewidth=1, label="Original: Average Difference")
    plt.plot(mod_maxs, linewidth=2, label="Modified: Max Difference")
    plt.plot(mod_means, linewidth=1, label="Modified: Mean Difference", color='orange')
    plt.xlabel("Time Units"); plt.ylabel("Difference")
    plt.xticks(np.arange(len(iters)), iters)
    plt.legend(loc="upper left")

    plt.subplot(122)
    plt.plot(org_norms, linewidth=3, label="Original: 2 Norm of Difference")
    plt.plot(mod_norms, linewidth=3, label="Modified: 2 Norm of Difference", color="red")
    plt.xlabel("Time Units"); plt.ylabel("Difference")
    plt.xticks(np.arange(len(iters)), iters)
    plt.legend(loc="upper left")

    plt.show()

analyze_difference(initial_diff1, initial_diff2, (500, 1000, 1500, 2000, 2500, 3000))
Original System --------------------------------------------------
After 500 units:
                Maximum difference:   0.2336170666411895
                Average difference:   0.02013832322734838
                2-Norm of difference: 0.8151417790156377
            
After 1000 units:
                Maximum difference:   22.164437248682955
                Average difference:   1.1682246287802942
                2-Norm of difference: 104.90295474163742
            
After 1500 units:
                Maximum difference:   25.433751805472287
                Average difference:   3.047811713432532
                2-Norm of difference: 234.30749671459355
            
After 2000 units:
                Maximum difference:   25.540634478290496
                Average difference:   4.225264422881091
                2-Norm of difference: 328.9242065287722
            
After 2500 units:
                Maximum difference:   25.540634478290496
                Average difference:   4.8741200365781205
                2-Norm of difference: 393.8593082924756
            
After 3000 units:
                Maximum difference:   25.540634478290496
                Average difference:   5.355396440083211
                2-Norm of difference: 449.4886075701734
            
Modified System --------------------------------------------------
After 500 units:
                Maximum difference:   0.47737145276164383
                Average difference:   0.04183464920953739
                2-Norm of difference: 1.7446706388482525
            
After 1000 units:
                Maximum difference:   9.522904066347227
                Average difference:   0.5066821488195269
                2-Norm of difference: 44.3299437863608
            
After 1500 units:
                Maximum difference:   16.275824085500556
                Average difference:   2.061942791352202
                2-Norm of difference: 153.7035971121914
            
After 2000 units:
                Maximum difference:   16.275824085500556
                Average difference:   2.4686997677714206
                2-Norm of difference: 187.12535239251233
            
After 2500 units:
                Maximum difference:   16.275824085500556
                Average difference:   3.047949499399024
                2-Norm of difference: 239.77659145837796
            
After 3000 units:
                Maximum difference:   16.360309307340053
                Average difference:   3.283827892951549
                2-Norm of difference: 275.8090389665492
            

Thus for our particular set of initial conditions, the modified system appears to be somewhat more stable than the original system with respect to perturbations in the initial conditions.

2.2: Varying the Solver Tolerances

We can also solve the system from the same initial point, but with different solver tolerances. To see any differences, we need to let the system evolve for a little longer than before.

In [7]:
# Solve the original and modified systems, with different solver tolerances.
tol_org1 = solve_original(x0=initial[:3], atol=1e-14, rtol=1e-12, t=50)
tol_org2 = solve_original(x0=initial[:3], atol=1e-15, rtol=1e-13, t=50)
tol_mod1 = solve_modified(x0=initial, atol=1e-14, rtol=1e-12, t=50)
tol_mod2 = solve_modified(x0=initial, atol=1e-15, rtol=1e-13, t=50)

tol_diff1 = np.abs(np.mean(tol_org1-tol_org2, axis=0))
tol_diff2 = np.abs(np.mean(tol_mod1-tol_mod2, axis=0))

# Plot results.
plot_double(tol_org1, tol_org2, color1='b', color2='g', title='Original System with Different Solver Tolerances')
plot_double(tol_mod1, tol_mod2, color1='r', color2='g', title='Modified System with Different Solver Tolerances')
In [8]:
analyze_difference(tol_diff1, tol_diff2, (1000, 2000, 3000, 4000, 5000))
Original System --------------------------------------------------
After 1000 units:
                Maximum difference:   9.063475863158033e-06
                Average difference:   4.101679623159729e-07
                2-Norm of difference: 4.367316496921059e-05
            
After 2000 units:
                Maximum difference:   0.06782569672457732
                Average difference:   0.0016954098542874504
                2-Norm of difference: 0.2755582003304529
            
After 3000 units:
                Maximum difference:   23.02726883812709
                Average difference:   1.4195116080066257
                2-Norm of difference: 219.11287802801306
            
After 4000 units:
                Maximum difference:   23.15718243077268
                Average difference:   2.955076135375699
                2-Norm of difference: 359.6158222125011
            
After 5000 units:
                Maximum difference:   23.15718243077268
                Average difference:   3.7150029262006488
                2-Norm of difference: 448.8535331636321
            
Modified System --------------------------------------------------
After 1000 units:
                Maximum difference:   7.18036890123841e-07
                Average difference:   1.5426836126466136e-08
                2-Norm of difference: 2.231996026510421e-06
            
After 2000 units:
                Maximum difference:   0.0003080385466248892
                Average difference:   9.636192833738549e-06
                2-Norm of difference: 0.0015986291848828092
            
After 3000 units:
                Maximum difference:   3.800343432216277
                Average difference:   0.06470419472398026
                2-Norm of difference: 15.987530084722696
            
After 4000 units:
                Maximum difference:   16.43239355344914
                Average difference:   1.0238754225365554
                2-Norm of difference: 177.2273683321611
            
After 5000 units:
                Maximum difference:   16.43239355344914
                Average difference:   1.7793372356941788
                2-Norm of difference: 256.79181451514984
            

Once again, it appears that the modified system is a little more stable than the original system (at least, for these initial conditions). However, this result must be taken with a grain of salt--it probably has more to do with the instabilities of the ODE solver $\textbf{scipy.integrate.odeint()}$ than the instabilities of the systems.

2.3: Perturbing the System Parameters

The parameters $b=\frac{8}{3}, \sigma=10, r=28$ give rise to the standard butterfly-like manifold that we have seen in all prior solutions of both the original and modified Lorenz systems. Changing one of the parameters slightly should result in a change in the solutions, as with the other perturbations we have used so far. We choose to vary $r$ from $27.9$ to $28.1$.

In [9]:
# Solve the original and modified systems, with different solver tolerances.
param_org1 = solve_original(x0=initial[:3], r=27.9)
param_org2 = solve_original(x0=initial[:3], r=28.1)
param_mod1 = solve_modified(x0=initial, r=27.9)
param_mod2 = solve_modified(x0=initial, r=28.1)

param_diff1 = np.abs(np.mean(param_org1-param_org2, axis=0))
param_diff2 = np.abs(np.mean(param_mod1-param_mod2, axis=0))

# Plot results.
plot_double(param_org1, param_org2, color1='b', color2='g', title='Original System with Different Parameters')
plot_double(param_mod1, param_mod2, color1='r', color2='g', title='Modified System with Different Parameters')
In [10]:
analyze_difference(param_diff1, param_diff2, (500, 1000, 1500, 2000))
Original System --------------------------------------------------
After 500 units:
                Maximum difference:   21.100089192062466
                Average difference:   1.741436853442377
                2-Norm of difference: 73.70347108115313
            
After 1000 units:
                Maximum difference:   23.72012814715785
                Average difference:   4.004748065538088
                2-Norm of difference: 209.53656535171825
            
After 1500 units:
                Maximum difference:   23.72012814715785
                Average difference:   4.768645405589466
                2-Norm of difference: 270.4458689574532
            
After 2000 units:
                Maximum difference:   23.72012814715785
                Average difference:   5.465712332070162
                2-Norm of difference: 337.9335741486734
            
Modified System --------------------------------------------------
After 500 units:
                Maximum difference:   13.335635859708503
                Average difference:   3.4696690167543687
                2-Norm of difference: 117.5432532498658
            
After 1000 units:
                Maximum difference:   13.335635859708503
                Average difference:   3.7321274037609156
                2-Norm of difference: 163.15723174185115
            
After 1500 units:
                Maximum difference:   14.660516018691515
                Average difference:   3.35097492808628
                2-Norm of difference: 186.20049774220524
            
After 2000 units:
                Maximum difference:   15.594759283721393
                Average difference:   3.411146860687478
                2-Norm of difference: 216.78557324582496
            

Once again, it looks like the modified system does a little better than the original system.

Section 2.4: Accounting for Extra Dimensions

What can explain these differences? Part of it has to do with the way that the ODE solver scipy.integrate.odeint() works, as already mentioned. The main factor, however, is the difference in dimensions. In fact, upon closer inspection, we see how the last two dimensions (corresponding the variables $U$ and $V$) might be affecting the analysis given above.

In [11]:
def examine(diff):
    for dim, dimension_array in zip(['X','Y','Z','U','V'], diff):
        print("\t2 norm of {} dimension of solutions:".format(dim), norm(dimension_array))

print("Perturbed initial conditions")
examine(np.abs(initial_mod1 - initial_mod2))
print("Different solver tolerances")
examine(np.abs(tol_mod1 - tol_mod2))
print("Perturbed system parameters")
examine(np.abs(param_mod1 - param_mod2))
Perturbed initial conditions
	2 norm of X dimension of solutions: 516.5853419837487
	2 norm of Y dimension of solutions: 632.9656177233671
	2 norm of Z dimension of solutions: 492.6604783050511
	2 norm of U dimension of solutions: 121.51657120690525
	2 norm of V dimension of solutions: 44.29920552491583
Different solver tolerances
	2 norm of X dimension of solutions: 502.7777373934655
	2 norm of Y dimension of solutions: 612.327723090078
	2 norm of Z dimension of solutions: 453.133729641105
	2 norm of U dimension of solutions: 119.68263718226515
	2 norm of V dimension of solutions: 41.07939551862201
Perturbed system parameters
	2 norm of X dimension of solutions: 442.8205682659842
	2 norm of Y dimension of solutions: 540.5807740177662
	2 norm of Z dimension of solutions: 357.7079501466286
	2 norm of U dimension of solutions: 106.01543973879178
	2 norm of V dimension of solutions: 33.04908525156817

In every case, the 2-norm of the $U$ and $V$ dimensions of the solutions are significantly smaller than the norms of the $X$, $Y$, and $Z$ dimensions. Before, when we took the mean the array of absolute differences, the 4th and 5th dimensions skewed the total difference to the left. This is problematic; the fact that the final two dimensions have smaller norms is actually not very significant because of the nondimensional nature of the system, so those numbers should not factor much into our results.

We fix this by repeating some of our previous steps, but now only using the first three dimensions of the solutions to the modified system.

In [12]:
# Reset the difference vectors by only averaging the first three dimensions.
initial_diff2 = np.abs(np.mean((initial_mod1-initial_mod2)[:3], axis=0))
tol_diff2 = np.abs(np.mean((tol_mod1-tol_mod2)[:3], axis=0))
param_diff2 = np.abs(np.mean((param_mod1-param_mod2)[:3], axis=0))

# Recreate the plots from the previous three subsections.
plt.suptitle("Perturbations in Initial Conditions", fontsize=20)
analyze_difference(initial_diff1, initial_diff2, (500, 1000, 1500, 2000, 2500, 3000), False)
plt.suptitle("Different Solver Tolerances", fontsize=20)
analyze_difference(tol_diff1, tol_diff2, (1000, 2000, 3000, 4000, 5000), False)
plt.suptitle("Perturbed System Parameters", fontsize=20)
analyze_difference(param_diff1, param_diff2, (500, 1000, 1500, 2000), False)

These new plots provide a new perspective: the 'stability' of the modified system is not actually much better than that of the original system. In fact, it only seems to be generally more stable relative to changes in the solver tolerances, which doesn't tell us much about the nature of the system itself. The original system actaully seems to be a little better than the modified system with respect to perturbation in the initial conditions and system parameters.

From this analysis it is fair to conclude that the modified Lorenz system is about as stable as the original system. That is, they are not significantly different. The 4th and 5th dimensions of the modified system seem to be somewhat well-behaved, but in $\mathbb{R}^3$ it behaves about as well (or rather, about as poorly) as the original system.

Section 3: Conditioning

Recall that the relative condition number $\kappa$ of a differentiable mapping $f: X \rightarrow Y$ between normed vector spaces ($f$ is referred to as the "problem") is given by

$$\kappa(\bf{x}) = \frac{\|\bf{J}(\bf{x})\|}{\|\bf{f}(\bf{x})\| / \|\bf{x}\|}$$

(we will continue using the 2-norm for the remainder of this project).

Section 3.1: Evaluating the Derivatives

For our two Lorenz systems, there are a few different conditioning problems that we can consider. The first is the problem of actually calculating the derivatives $\dot{X}$, $\dot{Y}$, $\dot{Z}$, $\dot{U}$ and $\dot{V}$ at a given point $(x,y,z,u,v)$. For the original system, we only use the first three dimensions.

For this problem, the Jacobian of the original system is given by $$J_{org}(x,y,z) = \left( \begin{array}{ccc} -\sigma & \sigma & 0 \\ r-Z & -1 & -X \\ Y & X & -b \end{array} \right)$$

and the Jacobian of the modified system is given by $$J_{mod}(x,y,z, u, v) = \left( \begin{array}{ccccc} -\sigma & \sigma & 0 & -p & 0\\ r-Z & -1 & -X & 0 & 0\\ Y & X & -b & 0 & 0\\ p-V & 0 & 0 & -\sigma & -X\\ U & 0 & 0 & X & -q\end{array} \right)$$

We will 'sample' the condition numbers at several points to try to estimate an upper bound on the condition number. In this case, we're actually estimating

$$\underset{\bf{x}\in\mathbb{R}^3}{\text{sup }}\kappa(\bf{x})$$
In [13]:
p = (np.sqrt(10)*(8/3.)**(1.5))/(.05*np.pi**2*2**(11./4.))

# Define functions for evaluating the condition number of this problem.
def val_org(x, y, z):
    """Calculate the value of the original Lorenz system at a point (with standard system parameters)."""
    return np.array([10*(y-x), 28*x - x*z - y, x*y - 8*z/3.])

def val_mod(x, y, z, u, v):
    """Calculate the value of the modified Lorenz system at a point (with standard system parameters)."""
    return np.array([10*(y-x) - p*u, 28*x - x*z - y, x*y - 8*z/3., -x*v + p*x - 10*u, x*u - 40*v/3.])

def Jac_val_org(x, y, z):
    """Return the Jacobian of the original Lorenz system at a point (with standard system parameters)."""
    return np.array([[ -10,  10,     0],
                     [28-z,  -1,     x],
                     [   y,   x, -8/3.]])

def Jac_val_mod(x, y, z, u, v):
    """Return the Jacobian of the modified Lorenz system at a point (with standard system parameters)."""
    return np.array([[ -10,  10,     0,  -p,      0],
                     [28-z,  -1,    -x,   0,      0],
                     [   y,   x, -8/3.,   0,      0],
                     [ p-v,   0,     0, -10,     -x],
                     [   u,   0,     0,   x, -40/3.]])

def cond_val_org(x,y,z):
    """Calculate the condition number of calculating the derivatives of x, y, and z at the given point."""
    return norm(Jac_val_org(x,y,z))/(norm(val_org(x,y,z))/norm(np.array([x,y,z])))
    
def cond_val_mod(x, y, z, u, v):
    """Calculate the condition number of calculating the derivatives of x, y, and z at the given point."""
    return norm(Jac_val_mod(x,y,z,u,v))/(norm(val_mod(x,y,z,u,v))/norm(np.array([x,y,z,u,v])))

# Estimate the condition number with 50,000 random points.
max_org_cond, max_mod_cond = 0, 0
for i in range(50000):
    point = np.random.random(5)*np.random.randint(-50,50,5)
    try:
        c1 = cond_val_org(point[0], point[1], point[2])
        c2 = cond_val_mod(point[0], point[1], point[2], point[3], point[4])
    except ZeroDivisionError:
        continue
    max_org_cond = max(max_org_cond, c1)
    max_mod_cond = max(max_mod_cond, c2)
print("Relative Condition Number for Original System: {}".format(max_org_cond),
      "Relative Condition Number for Modified System: {}".format(max_mod_cond), sep='\n')
Relative Condition Number for Original System: 60.875443314546594
Relative Condition Number for Modified System: 26.875438937389912

For this problem, it appears that the problem of calculating the derivatives in the modified system is actually a little more well-conditioned than the problem of calculating the derivatives in the original system. However, both condition numbers are pretty low, suggesting that the ODE solver does not perform significantly worse than usual. If the numbers were very high, we would have to seriously question the analysis from the previous section.

Section 3.2: Calculating the Fixed Points

To conclude, we estimate the condition number of another problem: that of finding the fixed points of the system.

Recall that a fixed point of a system $\bf{x}' = \bf{f}(\bf{x})$ is a point $\bf{\bar{x}}$ where all of the derivatives vanish. That is, we seek $\bf{\bar{x}}$ such that $\bf{f}(\bf{\bar{x}}) = \bf{0}$. Two of the three fixed points of both the original and modified systems in $\mathbb{R}^3$ are easy seen in the solution plots that we have seen so far as the attracting point on the solution manifold.

Section 3.2.1: Algebraic Solutions

In our case, finding the fixed points amounts to solving the nonlinear systems

$$\begin{align} \nonumber &0 = \sigma \left(Y-X\right) && 0 = \sigma \left(Y-X\right) - pU\\ \nonumber &0 = rX-XZ-Y && 0 = rX-XZ-Y\\ \nonumber &0 = XY-bZ && 0 = XY-bZ\\ \nonumber & && 0 = -XV + pX - \sigma U\\ \nonumber & && 0 = XU - qV \end{align}$$

Which can be reformulated as

$$\left( \begin{array}{ccc} -\sigma & \sigma & 0 \\ r & -1 & 0 \\ 0 & 0 & -b \end{array} \right) \left( \begin{array}{c}X\\Y\\Z\end{array}\right) + \left( \begin{array}{c}0\\-XZ\\XY\end{array}\right) = \left( \begin{array}{c}0\\0\\0\end{array}\right)$$

for the original system and

$$\left( \begin{array}{ccccc} -\sigma & \sigma & 0 & -p & 0\\ r & -1 & 0 & 0 & 0\\ 0 & 0 & -b & 0 & 0\\ p & 0 & 0 & -\sigma & 0\\ 0 & 0 & 0 & 0 & -q\end{array} \right) \left( \begin{array}{c}X\\Y\\Z\\U\\V\end{array}\right) + \left( \begin{array}{c}0\\-XZ\\XY\\-XV\\XU\end{array}\right) = \left( \begin{array}{c}0\\0\\0\end{array}\right)$$

for the modified system.

This is an awful problem to solve numerically! However, it is possible to get an explicit algebraic formula for the fixed points of each system. For the original system, we have the following:

$$0 = \sigma(Y - X) \Longrightarrow X = Y\\ 0 = XY-bZ \Longrightarrow Z = \frac{1}{b}X^2\\ 0 = rX-XZ-Y \Longrightarrow 0 = rX - \frac{1}{b}X^3 - X \Longrightarrow X = 0, \pm\sqrt{b(r-1)}$$

Then the three fixed points for the original system are $(0,0,0)$ and $(\pm\sqrt{b(r-1)}, \pm\sqrt{b(r-1)}, r-1)$, as long as $r>1$.

The explicit formula for the fixed points of the modified system is not nearly as clean. We do not show the complete derivation here, but a process similar to the one given above yields the following as the fixed points of the modified system:

$$X = \pm\sqrt{\frac{-\frac{p^2}{b} + \sigma\left(\frac{r-1}{q} - \frac{\sigma}{b}\right) \mp \sqrt{\left(\frac{p^2}{b} - \sigma(\frac{r-1}{q} - \frac{\sigma}{b})\right)^2 + \frac{4\sigma}{bq}\left(\sigma^2(r-1)-p^2\right)}}{\frac{2\sigma}{bq}}}\\ \\ \\ Y = \frac{brX}{X^2 + b}\\ Z = \frac{rX^2}{X^2 + b}\\ U = \frac{pqX}{X^2 + \sigma q}\\ V = \frac{pX^2}{X^2 + \sigma q}$$

Where $p = \frac{\sqrt{\sigma}b^{3/2}}{\epsilon\pi^2 2^{11/4}}$ and $q = (4-b)\sigma$ as usual.

Section 3.2.2: Problem Definition and Jacobian

For this problem, the inputs are the system parameters $b$, and $r$ for the original system (since the fixed points do not depend at all on $\sigma$), plus $\sigma$ and $\epsilon$ in the modified system. We restrict the problem to computing the fixed point in the octant $\{(x,y,z) : x\ge0,y\ge0,z\le0 \}$. For the original system, this is the function

$$f_{org}(b, r) = \left(\begin{array}{c}\sqrt{b(r-1)}\\ \sqrt{b(r-1)}\\r-1\end{array}\right)$$

Which has Jacobian $$J_{org}(b, r) = \left(\begin{array}{cc} \frac{r-1}{2\sqrt{b(r-1)}} & \frac{b}{2\sqrt{b(r-1)}}\\ \frac{r-1}{2\sqrt{b(r-1)}} & \frac{b}{2\sqrt{b(r-1)}}\\ 0 & 1\end{array}\right)$$

For the modified system, however, this is a much uglier problem. The function is given by

$$f_{mod}(b, r, \sigma, \epsilon) = \left(\begin{array}{c}X\\ \frac{brX}{X^2 + b}\\ \frac{rX^2}{X^2 + b}\\ \frac{pqX}{X^2 + \sigma q}\\ \frac{pX^2}{X^2 + \sigma q}\end{array}\right)$$

where

$$X = \sqrt{\frac{\sigma\left(\frac{r-1}{q} - \frac{\sigma}{b}\right) -\frac{p^2}{b} - \sqrt{\left(\sigma(\frac{r-1}{q} - \frac{\sigma}{b}) -\frac{p^2}{b}\right)^2 + \frac{4\sigma}{bq}\left(\sigma^2(r-1)-p^2\right)}}{\frac{2\sigma}{bq}}}$$

The Jacobian is quite difficult to construct analytically. Instead, we use the symmetric difference quotient as an esitmate for the derivatives:

$$f'(x) = \underset{h\rightarrow 0}{\text{ lim }}\frac{f(x+h) - f(x-h)}{2h}$$

which yields a good approximation $$f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}\text{ with }h = \sqrt{\epsilon_{machine}x}$$

In [14]:
def fix_org(b,r):
    """Calculate a fixed point of the original Lorenz system."""
    term = np.sqrt(b*(r-1))
    return np.array([term, term, r-1])

def fix_mod(b,r,sig,eps):
    """Calculate a fixed point of the modified Lorenz system."""
    p = (np.sqrt(sig)*(b)**(1.5))/(eps*np.pi**2*2**(11./4.))
    q = (4-b)*r
    
    A = sig*((r-1)/q - sig/b) - p**2/b
    B = ((4*sig)/(b*q))*((r-1)*sig**2 - p**2)
    
    X = (A - np.sqrt(A**2 + B))/((2*sig)/(b*q))
    
    C = r*X/(X**2 + b)
    D = p*X/(X**2 + sig*q)
    
    Y = b*C
    Z = X*C
    U = q*D
    V = X*D
    return np.array([X,Y,Z,U,V])

def Jac_fix_org(b,r):
    """The Jacobian for calculating a fixed point of the original Lorenz system."""
    numer = r-1
    denom = 2*np.sqrt(b*(numer))
    left = numer/denom
    right = b/denom
    return np.array([[left, right],[left, right],[0,1]])

def Jac_fix_mod(b,r,sig,eps):
    """Numerically estimate the Jacobian for calculating a fixed point of the modified Lorenz system."""
    
    # Define function for calculating X, Y, Z, U, and V.
    def P(b_,sig_,eps_):
        return (np.sqrt(sig_)*b_**(1.5))/(eps_*np.pi**2*2**(11./4.))
    def Q(b_,r_):
        return (4-b_)*r_
    def X(b_,r_,sig_,eps_):
        p = P(b_,sig_,eps_)
        q = Q(b_,r_)

        A = sig*((r-1)/q - sig/b) - p**2/b
        B = ((4*sig)/(b*q))*((r-1)*sig**2 - p**2)

        return (A - np.sqrt(A**2 + B))/((2*sig)/(b*q))
    
    def Y(b_,r_,sig_,eps_):
        x = X(b_,r_,sig_,eps_)
        return b_*r_*x/(x**2 + b_)
    
    def Z(b_,r_,sig_,eps_):
        x = X(b_,r_,sig_,eps_)
        return r_*x**2/(x**2 + b)
    
    def U(b_,r_,sig_,eps_):
        p = P(b_,sig_,eps_)
        q = Q(b_,r_)
        x = X(b_,r_,sig_,eps_)
        return p*q*x/(x**2 + sig_*q)

    def V(b_,r_,sig_,eps_):
        p = P(b_,sig_,eps_)
        q = Q(b_,r_)
        x = X(b_,r_,sig_,eps_)
        return p*x**2/(x**2 + sig_*q)

    # Estimate the Jacobian using the symmetric difference quotient.
    h_b = 1e-8*np.sqrt(b)
    h_r = 1e-8*np.sqrt(r)
    h_s = 1e-8*np.sqrt(sig)
    h_e = 1e-8*np.sqrt(eps)

    dXdb = (X(b+h_b,r,sig,eps) - X(b-h_b,r,sig,eps))/(2*h_b)
    dXdr = (X(b,r+h_r,sig,eps) - X(b,r-h_r,sig,eps))/(2*h_r)
    dXds = (X(b,r,sig+h_s,eps) - X(b,r,sig-h_s,eps))/(2*h_s)
    dXde = (X(b,r,sig,eps+h_e) - X(b,r,sig,eps-h_e))/(2*h_e)
    
    dYdb = (Y(b+h_b,r,sig,eps) - Y(b-h_b,r,sig,eps))/(2*h_b)
    dYdr = (Y(b,r+h_r,sig,eps) - Y(b,r-h_r,sig,eps))/(2*h_r)
    dYds = (Y(b,r,sig+h_s,eps) - Y(b,r,sig-h_s,eps))/(2*h_s)
    dYde = (Y(b,r,sig,eps+h_e) - Y(b,r,sig,eps-h_e))/(2*h_e)
    
    dZdb = (Z(b+h_b,r,sig,eps) - Z(b-h_b,r,sig,eps))/(2*h_b)
    dZdr = (Z(b,r+h_r,sig,eps) - Z(b,r-h_r,sig,eps))/(2*h_r)
    dZds = (Z(b,r,sig+h_s,eps) - Z(b,r,sig-h_s,eps))/(2*h_s)
    dZde = (Z(b,r,sig,eps+h_e) - Z(b,r,sig,eps-h_e))/(2*h_e)
    
    dUdb = (U(b+h_b,r,sig,eps) - U(b-h_b,r,sig,eps))/(2*h_b)
    dUdr = (U(b,r+h_r,sig,eps) - U(b,r-h_r,sig,eps))/(2*h_r)
    dUds = (U(b,r,sig+h_s,eps) - U(b,r,sig-h_s,eps))/(2*h_s)
    dUde = (U(b,r,sig,eps+h_e) - U(b,r,sig,eps-h_e))/(2*h_e)
    
    dVdb = (V(b+h_b,r,sig,eps) - V(b-h_b,r,sig,eps))/(2*h_b)
    dVdr = (V(b,r+h_r,sig,eps) - V(b,r-h_r,sig,eps))/(2*h_r)
    dVds = (V(b,r,sig+h_s,eps) - V(b,r,sig-h_s,eps))/(2*h_s)
    dVde = (V(b,r,sig,eps+h_e) - V(b,r,sig,eps-h_e))/(2*h_e)

    return np.array([[dXdb, dXdr, dXds, dXde],
                     [dYdb, dYdr, dYds, dYde],
                     [dZdb, dZdr, dZds, dZde],
                     [dUdb, dUdr, dUds, dUde],
                     [dVdb, dVdr, dVds, dVde]])
    

def cond_fix_org(b,r):
    """Calculate the condition number for calculating a fixed point of the original Lorenz system."""
    return norm(Jac_fix_org(b,r))/(norm(fix_org(b,r))/norm(np.array([b,r])))

def cond_fix_mod(b,r,sig,eps):
    """Calculate the condition number for calculating a fixed point of the original Lorenz system."""
    return norm(Jac_fix_mod(b,r,sig,eps))/(norm(fix_mod(b,r,sig,eps))/norm(np.array([b,r,sig,eps])))


# Calculate the condition number with standard system parameters.
print("Relative Condition Number for Original System:",cond_fix_org(8/3.,28))
print("Relative Condition Number for Modified System:",cond_fix_mod(8/3.,28,10,.05))
Relative Condition Number for Original System: 2.35341829974
Relative Condition Number for Modified System: 156.30532443

For the original system, the condition number is fantastic. For the modified system, the condition number isn't terrible, but it's much worse than that of the original system. What happened?

We may expect that the larger number of dimensions in the modified system will inflate the condition number somewhat. Also, because the Jacobian is only an approximation (and probably a pretty poor one), we may take the condition number with a grain of salt. However, it is clear that the problem of calculating the fixed points of the modified system is a little more ill-conditioned than the problem of calculating the fixed points of the original system.

Section 4: Conclusion

The general result we can pull from these two analyses is that our modified Lorenz system is not significantly worse (or better) than the original Lorenz system in terms of stability and conditioning. This may seem like a boring result, but it is actually good to know that the modified system behaves similarly to the original system, indicating that the modified system is a legitimate variant of the original.